## estimates an individual-level regession of 2016 voting against 2012 voting plus controls,
## for voters in counties where either no cuts were made, or where hours were reduced

library(foreign)
library(lmtest)
library(effects)
library(nnet)
library(MASS)
library(tidyverse)
library(plyr) 
library(dplyr)
library(haven)
library(simcf)
library(ggplot2)
library(multiwayvcov)
library(stargazer)
options(scipen=20)



## attach the vbf object unless it is already attached
if (!length(grep("vbf_output", search()))) {
  cat ("Attaching vbf... this may take a few moments...")
  attach("vbf_output")
}

vbf$race <- as.factor(vbf$race)
vbf$age_category <- as.factor(vbf$age_category)
vbf <- within(vbf, race <- relevel(race, ref = "white"), age_category <- relevel(age_category, ref = "age_middle"))

# Function to calculate clustered SEs 
robustse.f <- function(model, cluster, df_correction) {
  require(sandwich)
  require(lmtest)
  require(multiwayvcov)
  if(missing(cluster)) {
    name <- deparse(substitute(model))
    modelname <- paste(name,"rob",sep=".")
    model$se <- coeftest(model, vcov=vcovHC(model,"HC1"))[,2]
    model$vcovHC <- vcovHC(model,"HC1")
    assign(modelname,model,envir = .GlobalEnv)
    coeftest(model, vcov=vcovHC(model,"HC1"))
  } else {
    name <- deparse(substitute(model))
    modelname <- paste(name,"rob",sep=".")
    vcovCL <- cluster.vcov(model, cluster, df_correction = df_correction)
    model$vcovCL <- vcovCL
    modelname <- paste(name,"clustrob",sep=".")
    model$se <- coeftest(model, vcovCL)[,2]
    assign(modelname,model,envir = .GlobalEnv)
    coeftest(model, vcovCL)
  }
}

# identify the counties that either reduced TOTAL HOURS or made no changes
total_hours_reduced<-data.frame(filter(vbf, totalhoursChg <=0))

summary(model_allchanges <- glm(factor(outcome2016 != "novote2016") ~ totalhoursChg + evehoursChg + 
                                  sathoursChg+sunhoursChg+sitesChg+
                                  factor(race) + factor(female) + factor(party) + 
                                  factor(age_category) + factor(outcome2012 != "novote2012"),
                                 data = total_hours_reduced %>% filter (voter_status_desc_2012 == "ACTIVE",
                                                       voter_status_desc_2016 == "ACTIVE",
                                                       ethnicity == "latino_not_yes",
                                                       (race == "black" | race == "white"),
                                                       (party == "dem" | party=="rep")),
                                family = "binomial"))


model_totalhoursChg_cluster <- robustse.f(model_allchanges, model_allchanges$data$COUNTY_NAM, df_correction = TRUE)

# identify the counties that either reduced EVENING HOURS or made no changes
evening_hours_reduced<-data.frame(filter(vbf, evehoursChg <=0))

summary(model_allchanges <- glm(factor(outcome2016 != "novote2016") ~ totalhoursChg + evehoursChg + 
                                  sathoursChg+sunhoursChg+sitesChg+
                                  factor(race) + factor(female) + factor(party) + 
                                  factor(age_category) + factor(outcome2012 != "novote2012"),
                                data = evening_hours_reduced %>% filter (voter_status_desc_2012 == "ACTIVE",
                                                                       voter_status_desc_2016 == "ACTIVE",
                                                                       ethnicity == "latino_not_yes",
                                                                       (race == "black" | race == "white"),
                                                                       (party == "dem" | party=="rep")),
                                family = "binomial"))


model_evehoursChg_cluster <- robustse.f(model_allchanges, model_allchanges$data$COUNTY_NAM, df_correction = TRUE)

# identify the counties that either reduced SATURDAY HOURS or made no changes
saturday_hours_reduced<-data.frame(filter(vbf, sathoursChg <=0))

summary(model_allchanges <- glm(factor(outcome2016 != "novote2016") ~ totalhoursChg + evehoursChg + 
                                  sathoursChg+sunhoursChg+sitesChg+
                                  factor(race) + factor(female) + factor(party) + 
                                  factor(age_category) + factor(outcome2012 != "novote2012"),
                                data = saturday_hours_reduced %>% filter (voter_status_desc_2012 == "ACTIVE",
                                                                         voter_status_desc_2016 == "ACTIVE",
                                                                         ethnicity == "latino_not_yes",
                                                                         (race == "black" | race == "white"),
                                                                         (party == "dem" | party=="rep")),
                                family = "binomial"))


model_sathoursChg_cluster <- robustse.f(model_allchanges, model_allchanges$data$COUNTY_NAM, df_correction = TRUE)

# identify the counties that either reduced SUNDAY HOURS or made no changes
sunday_hours_reduced<-data.frame(filter(vbf, sunhoursChg <=0))

summary(model_allchanges <- glm(factor(outcome2016 != "novote2016") ~ totalhoursChg + evehoursChg + 
                                  sathoursChg+sunhoursChg+sitesChg+
                                  factor(race) + factor(female) + factor(party) + 
                                  factor(age_category) + factor(outcome2012 != "novote2012"),
                                data = sunday_hours_reduced %>% filter (voter_status_desc_2012 == "ACTIVE",
                                                                          voter_status_desc_2016 == "ACTIVE",
                                                                          ethnicity == "latino_not_yes",
                                                                          (race == "black" | race == "white"),
                                                                          (party == "dem" | party=="rep")),
                                family = "binomial"))


model_sunhoursChg_cluster <- robustse.f(model_allchanges, model_allchanges$data$COUNTY_NAM, df_correction = TRUE)

# identify the counties that either reduced SITES or made no changes
sites_reduced<-data.frame(filter(vbf, sitesChg <=0))

summary(model_allchanges <- glm(factor(outcome2016 != "novote2016") ~ totalhoursChg + evehoursChg + 
                                  sathoursChg+sunhoursChg+sitesChg+
                                  factor(race) + factor(female) + factor(party) + 
                                  factor(age_category) + factor(outcome2012 != "novote2012"),
                                data = sites_reduced %>% filter (voter_status_desc_2012 == "ACTIVE",
                                                                        voter_status_desc_2016 == "ACTIVE",
                                                                        ethnicity == "latino_not_yes",
                                                                        (race == "black" | race == "white"),
                                                                        (party == "dem" | party=="rep")),
                                family = "binomial"))


model_sitesChg_cluster <- robustse.f(model_allchanges, model_allchanges$data$COUNTY_NAM, df_correction = TRUE)


# stargaze and write to tex file

stargazer(model_totalhoursChg_cluster, model_evehoursChg_cluster,model_sathoursChg_cluster,
          model_sunhoursChg_cluster,model_sitesChg_cluster,
          no.space = TRUE, 
          font.size="footnotesize",
          title = "Reductions in EIP availability and 2016 voter turnout",
          dep.var.labels.include = FALSE,
          covariate.labels=c("Percent Change: Total Hours",
                             "Percent Change: Evening Hours",
                             "Percent Change: Saturday Hours",
                             "Percent Change: Sunday Hours",
                             "Percent Change: EIP Sites",
                             "Black","Female","Republican",
                             "Age: 18-34","Age: 65+","Voted in 2012",
                             "Constant"),
          order=c(1,2,3,4,5,6,7,8,10,9,11,12),
          label = "tab:reduced_base_regression",
          model.numbers = FALSE,
          column.labels = c("Total Hours",
                            "Evening Hours",
                            "Saturday Hours",
                            "Sunday Hours",
                            "Sites"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "models_Chg_cluster_reduced.tex",
          dep.var.caption = "")


###################################################
######### INTERACTION WITH PARTY AND RACE #########
###################################################

# TOTAL HOURS

summary(model_allchanges <- glm(factor(outcome2016 != "novote2016") ~ totalhoursChg*factor(race)*factor(party) + evehoursChg + 
                                  sathoursChg+sunhoursChg+sitesChg+
                                   factor(female) + factor(age_category) + factor(outcome2012 != "novote2012"),
                                data = total_hours_reduced %>% filter (voter_status_desc_2012 == "ACTIVE",
                                                                       voter_status_desc_2016 == "ACTIVE",
                                                                       ethnicity == "latino_not_yes",
                                                                       (race == "black" | race == "white"),
                                                                       (party == "dem" | party=="rep")),
                                family = "binomial"))


model_totalhoursChg_cluster <- robustse.f(model_allchanges, model_allchanges$data$COUNTY_NAM, df_correction = TRUE)

stargazer(model_totalhoursChg_cluster,
          no.space = TRUE, 
          title = "Reductions in total EIP hours, race, party and 2016 voter turnout",
          dep.var.labels.include = TRUE,
          order=c(12,13,14,15,1,2,3,4,5,6,7,8,9,10,11),
          covariate.labels=c("Total Hours*Black","Total Hours*Republican",
                             "Black*Republican","Total Hours*Black*Republican",
                             "Total Hours","Black","Republican",
                             "Evening Hours","Saturday Hours","Sunday Hours","Sites",
                             "Female","Age: 18-34","Age: 65+","Voted 2012","Constant"),
          label = "tab:reduced_total_interact_regression",
          model.numbers = FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = c("Reduced Total Hours"),
          out = "totalhoursChg_cluster_interact_reduced.tex",
          dep.var.caption = "")

# EVENING HOURS 

summary(model_allchanges <- glm(factor(outcome2016 != "novote2016") ~ totalhoursChg + evehoursChg*factor(race)*(party) + 
                                  sathoursChg+sunhoursChg+sitesChg+factor(female) + 
                                  factor(age_category) + factor(outcome2012 != "novote2012"),
                                data = evening_hours_reduced %>% filter (voter_status_desc_2012 == "ACTIVE",
                                                                         voter_status_desc_2016 == "ACTIVE",
                                                                         ethnicity == "latino_not_yes",
                                                                         (race == "black" | race == "white"),
                                                                         (party == "dem" | party=="rep")),
                                family = "binomial"))


model_evehoursChg_cluster <- robustse.f(model_allchanges, model_allchanges$data$COUNTY_NAM, df_correction = TRUE)

stargazer(model_evehoursChg_cluster,
          no.space = TRUE, 
          title = "Reductions in evening EIP hours, race, party and 2016 voter turnout",
          dep.var.labels.include = TRUE,
          order=c(12,13,14,15,1,2,3,4,5,6,7,8,9,10,11),
          covariate.labels=c("Evening Hours*Black","Evening Hours*Republican",
                             "Black*Republican","Evening Hours*Black*Republican",
                             "Evening Hours","Black","Republican",
                             "Total Hours","Saturday Hours","Sunday Hours","Sites",
                             "Female","Age: 18-34","Age: 65+","Voted 2012","Constant"),
          label = "tab:reduced_eve_interact_regression",
          model.numbers = FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = c("Reduced Evening Hours"),
          out = "evehoursChg_cluster_interact_reduced.tex",
          dep.var.caption = "")

# SATURDAY HOURS 

summary(model_allchanges <- glm(factor(outcome2016 != "novote2016") ~ totalhoursChg + evehoursChg + 
                                  sathoursChg*factor(race)*factor(party)+
                                  sunhoursChg+sitesChg+factor(female) + 
                                  factor(age_category) + factor(outcome2012 != "novote2012"),
                                data = saturday_hours_reduced %>% filter (voter_status_desc_2012 == "ACTIVE",
                                                                          voter_status_desc_2016 == "ACTIVE",
                                                                          ethnicity == "latino_not_yes",
                                                                          (race == "black" | race == "white"),
                                                                          (party == "dem" | party=="rep")),
                                family = "binomial"))


model_sathoursChg_cluster <- robustse.f(model_allchanges, model_allchanges$data$COUNTY_NAM, df_correction = TRUE)

stargazer(model_sathoursChg_cluster,
          no.space = TRUE, 
          title = "Reductions in Saturday EIP hours, race, party and 2016 voter turnout",
          dep.var.labels.include = TRUE,
          order=c(12,13,14,15,1,2,3,4,5,6,7,8,9,10,11),
          covariate.labels=c("Saturday Hours*Black","Saturday Hours*Republican",
                             "Black*Republican","Saturday Hours*Black*Republican",
                             "Saturday Hours","Black","Republican",
                             "Total Hours","Evening Hours","Sunday Hours","Sites",
                             "Female", "Age: 18-34","Age: 65+","Voted 2012","Constant"),
          label = "tab:reduced_sat_interact_regression",
          model.numbers = FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = c("Reduced Saturday Hours"),
          out = "sathoursChg_cluster_interact_reduced.tex",
          dep.var.caption = "")


# SUNDAY HOURS 
summary(model_allchanges <- glm(factor(outcome2016 != "novote2016") ~ totalhoursChg + evehoursChg + 
                                  sathoursChg+sunhoursChg*factor(race)*factor(party)+
                                  sitesChg+factor(female) +factor(age_category) + factor(outcome2012 != "novote2012"),
                                data = sunday_hours_reduced %>% filter (voter_status_desc_2012 == "ACTIVE",
                                                                        voter_status_desc_2016 == "ACTIVE",
                                                                        ethnicity == "latino_not_yes",
                                                                        (race == "black" | race == "white"),
                                                                        (party == "dem" | party=="rep")),
                                family = "binomial"))


model_sunhoursChg_cluster <- robustse.f(model_allchanges, model_allchanges$data$COUNTY_NAM, df_correction = TRUE)

stargazer(model_sunhoursChg_cluster,
          no.space = TRUE, 
          title = "Reductions in Sunday EIP hours, race, party and 2016 voter turnout",
          dep.var.labels.include = TRUE,
          order=c(12,13,14,15,1,2,3,4,5,6,7,8,9,10,11),
          covariate.labels=c("Sunday Hours*Black","Sunday Hours*Republican",
                             "Black*Republican","Sunday Hours*Black*Republican",
                             "Sunday Hours","Black","Republican",
                             "Total Hours","Evening Hours","Saturday Hours","Sites",
                             "Female","Age: 18-34","Age: 65+","Voted 2012","Constant"),
          label = "tab:reduced_sun_interact_regression",
          model.numbers = FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = c("Reduced Sunday Hours"),
          out = "sunhoursChg_cluster_interact_reduced.tex",
          dep.var.caption = "")

# SITES 

summary(model_allchanges <- glm(factor(outcome2016 != "novote2016") ~ totalhoursChg + evehoursChg + 
                                  sathoursChg+sunhoursChg+sitesChg*factor(race)*factor(party)+
                                  factor(female) +factor(age_category) + factor(outcome2012 != "novote2012"),
                                data = sites_reduced %>% filter (voter_status_desc_2012 == "ACTIVE",
                                                                 voter_status_desc_2016 == "ACTIVE",
                                                                 ethnicity == "latino_not_yes",
                                                                 (race == "black" | race == "white"),
                                                                 (party == "dem" | party=="rep")),
                                family = "binomial"))

model_sitesChg_cluster <- robustse.f(model_allchanges, model_allchanges$data$COUNTY_NAM, df_correction = TRUE)

stargazer(model_sitesChg_cluster,
          no.space = TRUE, 
          title = "Changes in EIP sites, race, party and 2016 voter turnout",
          dep.var.labels.include = TRUE,
          order=c(12,13,14,15,1,2,3,4,5,6,7,8,9,10,11),
          covariate.labels=c("Sites Hours*Black","Sites Hours*Republican",
                             "Black*Republican","Sites Hours*Black*Republican",
                             "Sites Hours","Black","Republican",
                             "Total Hours","Evening Hours","Saturday Hours","Sunday Hours",
                             "Female", "Age: 18-34","Age: 65+","Voted 2012","Constant"),
          label = "tab:reduced_sites_interact_regression",
          model.numbers = FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = c("Reduced Sites"),
          out = "sitesChg_cluster_interact_reduced.tex",
          dep.var.caption = "")





